Modulational instability of two pairs of counter-propagating 
waves and energy exchange in two-component media 

S.D. Griffiths , R.H.J. Grimshaw 6 , K.R. Khusnutdinova 6 * 

a Department of Atmospheric Sciences, University of Washington, 
Box 351640, Seattle, WA 98195, USA 

department of Mathematical Sciences, Loughborough University, 
Loughborough, LE11 3TU, UK 

Abstract 

The dynamics of two pairs of counter-propagating waves in two- component media is considered 
within the framework of two generally nonintegrable coupled Sine-Gordon equations. We consider 
the dynamics of weakly nonlinear wave packets, and using an asymptotic multiple-scales expansion 
we obtain a suite of evolution equations to describe energy exchange between the two components 
of the system. Depending on the wave packet length-scale vis-a-vis the wave amplitude scale, these 
evolution equations are either four non-dispersive and nonlinearly coupled envelope equations, or 
four non-locally coupled nonlinear Schrodinger equations. We also consider a set of fully coupled 
nonlinear Schrodinger equations, even though this system contains small dispersive terms which 
are strictly beyond the leading order of the asymptotic multiple-scales expansion method. Using 
both the theoretical predictions following from these asymptotic models and numerical simulations 
of the original unapproximated equations, we investigate the stability of plane- wave solutions, and 
show that they may be modulationally unstable. These instabilities can then lead to the formation 
of localized structures, and to a modification of the energy exchange between the components. 
When the system is close to being integrable, the time- evolution is distinguished by a remarkable 
almost periodic sequence of energy exchange scenarios, with spatial patterns alternating between 
approximately uniform wavetrains and localized structures. 
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1 Introduction 

The dynamics of nonlinear waves is central to understanding the behaviour of a wide range of physical 
systems. Applications include Langmuir waves in plasmas (e.g. 1 ), waves in optical fibres (e.g. 
J2J 01) , and water waves (e.g. 4 and references therein). An important theme of much of this work 
has been understanding modulational instabilities of plane waves, and in recent years a particular 
focus has been the instability of nonlinear ly coupled waves. Examples include studies of counter- 
propagating waves in parametrically forced systems (e.g. HOEOH]), standing water waves (e.g. jSJ), 
co-propagating waves in shallow water j^j, and more general energy exchange processes in coupled 
Sine-Gordon equations [U3- Other interesting recent results related to the dynamics of nonlinearly 
coupled waves include the rigorous proof of the existence of standing water waves [ll] , and a study of 
resonant triad dynamics in weakly damped Faraday waves with two- frequency forcing • 

Our present study is devoted to multi-component systems that may be modelled by nonlinear 
partial differential equations with one space dimension x, and with time t, which admit linear stable 
wave solutions of the form Aexp i(kx — cot). Here uo = uo(k), where importantly this dispersion relation 
has several branches. We aim to study the nonlinear development of several plane waves with the same 
wavenumber /c, taking particular care to diagnose any modulational instabilities. As is customary we 
shall assume that each wave amplitude A is 0(e) where £ < 1, thus enabling the weakly nonlinear 
regime to be studied using asymptotic methods. 

When there is just a single wave, the procedure is well-known. There is a distinguished lengthscale 
for modulations to the wave amplitude, 0(e _1 ), on which the effects of dispersion and nonlinearity 
become simultaneously important, on a timescale 0(e~ 2 ). Then, in a reference frame moving with the 
group velocity of the wave, the governing nonlinear partial differential equations can be reduced via 
the method of multiple scales to a nonlinear Schrodinger (NLS) equation for the wave amplitude A 
(e.g. ^2j). As is well known, if uo" (k) has the same sign as the coefficient of the cubic nonlinear term, 
then a plane wave is modulationally unstable 14 (sometimes called a Benjamin-Feir instability, or 
sideband instability [Ej). 

One can apply a similar methodology for the case of two waves, possibly counter-propagating, say 
Aex.pi(kx — uoi(k)t) + B ex.pi(kx — uj2{k)t). With the same scaling as above, in general the governing 
nonlinear partial differential equations can be reduced to non-locally coupled NLS equations for A 
and B [HOElIIBj. For unbounded or periodic systems, the coupling terms either disappear or can 
be eliminated by a phase transformation 17 , and so modulational instability is determined by the 
instability of the individual waves. However, if the group velocities uj[ and uo' 2 are equal, or differ only 
by an 0(e) quantity, then one obtains locally coupled NLS equations, and there may be additional 
instabilities solely due to this coupling, even when both waves are individually stable pi 120) . 

More generally, one can pursue an asymptotic multiple-scales expansion in which the wave ampli- 
tudes, of 0(e), are modulated on a lengthscale of 0(a _1 ). Then, retaining all terms of 0(ea, ea 2 , £ 3 ), 
the governing partial differential equations can be reduced to coupled NLS equations for the two wave 
amplitudes. Expressed here in unsealed variables, they are of the form 

A t + v 91 A x = (iu'{/2)A xx + im\A\ 2 A + i^\B\ 2 A, (1) 

where v 91 = co[ is the group velocity, and \±\^ are constant nonlinear coefficients. There is a corre- 
sponding equation for £?, with a different transport term due to the group velocity v 92 = uo 2 . On 
the right-hand side, the linear dispersive terms of 0(ecr 2 ) are retained, even though these are strictly 
smaller than the O(ecr) transport terms on the left-hand side. Note that for a single wave, the trans- 
port term can be removed by moving to the group velocity reference frame, but this is generally not 
possible in coupled systems with differing group velocities. The nonlinear terms on the right-hand 
side are 0(e 3 ), and with the conventional NLS choice that a = e, these equations can then be system- 
atically reduced to the aforementioned non-locally coupled NLS equations, provided that v 91 is not 
equal or close to v 92 . Alternatively, one could choose a = e 2 giving a balance between the transport 
terms and the nonlinear terms, leading to non-dispersive and nonlinearly coupled envelope equations. 
Nevertheless, although coupled NLS equations of the form are not generally asymptotically exact, 
in that they contain small dispersive terms which are strictly beyond the leading order transport 
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terms, they may be useful as they do capture the leading order dispersive and nonlinear effects. Stud- 
ies based on this type of approach have recently been made for the parametric generation of waves 
by an external oscillatory field in the framework of two coupled, damped, parametrically driven NLS 
equations (see 01110 and references therein), and for the interaction of two waves in shallow water 
in the framework of the KdV equation 0. 

The aims of this study are to extend our understanding of the nonlinear behaviour of multiple- 
wave systems, to derive asymptotic equations to describe many interacting waves, and to assess the 
applicability and accuracy of the various resulting coupled wave models, with particular reference to 
the modulational instability of plane waves. This work is performed in the context of energy exchange 
processes in two-component systems, permitting us to examine explicitly both two- wave and four- wave 
dynamics. We consider a system of coupled nonlinear Klein-Gordon equations 

Utt ~ Uxx = f u (u, w), w tt - c 2 w xx f w (u, w), (2) 

where the subscripts denote partial derivatives. Here, f(u, w) is a potential function for the nonlinear 
coupling, and c is the ratio of the acoustic velocities of the components u and w. In its simplest form, 
0) describes the long-wave dynamics of two coupled one-dimensional periodic chains of particles, the 
elements of each chain being linked by linear coupling, and the chains themselves being linked by the 
nonlinear coupling [21] . Applications of this model include nonlinear waves in bi-layers with imperfect 
interfaces [22j [22) , and dynamical processes in DNA [2] HHj • 

The particular case we shall study is that for which the potential function is f(u,w) = cos(5u — 
w) — 1. Then the system (0 reduces to a set of coupled Sine-Gordon equations: 

u t t - u xx = -S 2 sin(w - w), w tt - c 2 w xx = sin(u - w), (3) 

where the variable u replaces 5u, compared to (0. In addition to the two applications mentioned 
above, the coupled Sine- Gordon equations generalize the Frenkel-Kontorova dislocation model 0J 
(see also 27 , and references therein). As shown in [2EJ, processes involving an exchange of energy 
between the two components of the system are possible due to the existence of two branches of 
the dispersion relation. These processes represent a continuum generalization of energy exchange in a 
system of coupled pendulums 29 (see also [SJ). The mechanism of such energy exchange is essentially 
linear, with energy being transferred between the physical components of the system, rather than the 
wave components. Therefore, a natural question to ask is, what happens to these linear plane wave 
solutions under the influence of nonlinearity, especially if they are modulationally perturbed? Here, 
we are interested in studying such modulational instabilities in a more or less generic situation, in 
contrast to the recent papers on targeted energy transfer, which placed certain restrictions on the 
physical properties of the system |21 E3 • 

First steps in this direction have been made in our recent work • There we dealt with a particular 
solution of (0 involving only a single pair of waves, and considered its modulational instability in the 
context of two non-locally coupled NLS equations. In the present study, we generalise this analysis to 
two pairs of counter-propagating waves, which generally leads to a system of four non-locally coupled 
NLS equations. Note that the four- wave case is essentially different from the two- wave case, since the 
number of characteristic variables now differs from the number of independent variables, making the 
strict derivation of the asymptotic evolution equations not so straightforward. However, as we show 
here, the same methodology can, in general, be used to derive n non-locally coupled NLS equations 
for the co-evolution of n waves. We also re-examine the two- wave stability problem using the coupled 
NLS equations, and compare the predictions with those of the non-locally coupled NLS equations. 
Throughout, numerical simulations of the governing coupled Sine- Gordon equations are used to assess 
the various predictions over a range of wave amplitudes. 

Our paper is organized as follows. In Section we discuss solutions exhibiting periodic energy 
exchange between the two components. Section is devoted to the derivation of asymptotic models 
for weakly nonlinear waves. Specifically, we derive non-locally coupled NLS equations for the regime 
with modulations on a lengthscale of 0(e _1 ), and also non-dispersive nonlinear ly coupled envelope 
equations for the regime with modulations on a lengthscale of 0(e~ 2 ). Then we motivate the intro- 
duction of the coupled NLS equations, and give them for the case of two pairs of counter-propagating 
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waves (correcting an error in 28 ). In Section QJ we study modulational instability of weakly nonlin- 
ear coupled two- wave solutions in each of the respective models, extending our previous study [TU] . 
In Section [SJ we study modulational instability of weakly nonlinear four- wave solutions, using two 
asymptotic models. Numerical simulations of the system (0 of coupled Sine- Gordon equations are 
presented in Section^ with a focus on the growth of any modulational instabilities, and the numerical 
results are compared with the theoretical predictions of Sections 0] and We conclude in Section [3 



2 Energy exchange 

The system (0 is Lagrangian, with density 

L = ^(u 2 + w 2 t -u\- c 2 w 2 x ) + f(u, w). 

For any function f(u, w), the system (0 has conservation laws for energy and momentum. (In partic- 
ular cases, there are additional conservation laws, which can be found using results from a Lie group 
classification of these equations [221, an d by using the Noether theorem (e.g. 
The energy conservation law for the system (0 has the form 

Q (u 2 t + u 2 x + S 2 w 2 + c 2 S 2 w 2 x ) + S 2 (1 - cos(u - w))^j - -|- (u t u x + c 2 S 2 w t w x ) = 0. 

For periodic solutions in < x < /, we define domain integrated average energies 

Eu = l J l + u *) dx ' Ew = 1 1 T + ° 2w ^ dx ' (4) 
i 

E c = - I 5 2 (l-cos(u-w))dx, (5) 



describing the energy in the ^-component, the energy in the ^-component, and the energy of coupling, 
respectively. Then the domain integrated average energy E is simply 

E = E U + E W + E C . (6) 

Although E is constant, we can identify solutions that have a time-periodic exchange of energy between 
E u , E w and E c . 

2.1 Linear solutions 

When \u — w\ <C 1, the coupled Sine-Gordon equations (0 reduce to a linear system. Looking for a 
solution of the form (u, w) = (1, a) ex.pi(kx — uot) + c.c. leads to the dispersion relation 



'1,2 - g 



(7) 



where v\ = 5 2 + & 2 and z/| = 1 + c 2 k 2 . The corresponding values of a, giving the ratio w/iz, are 



"1,2 



S<2 ^2-^1,2' 

We note that ai > 0, and ct2 < 0. For example, if c = 1, then a\ = 1 and «2 = — £ -2 . 

We first consider a superposition of two linear waves, taken to be right-propagating, with the 
same positive wavenumber k but different frequencies coi(k) and ^(fe). Further, we suppose that one 
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component, taken to be u, is initially excited, while the other component w is initially at rest. The 
corresponding initial conditions are 



u(x, 0) = U cos kx, uAx, 0) = ^ 1 \ U20Ll ^ jj sm fc x , 

Oil + \Oi2\ 

( n x n / n x (^l -^2)ai|«2| rr • 7 
im:r,(JJ = 0, i^(x,(J) = : : U sin/cx, 

<*1 + \OL2\ 



(8) 



and one may write the time-dependent solution as 

u = U(t) cos(kx — 7+£ — 0(t)), w = ai J7 shi7_£ sin(kx — 7+t), (9) 

CKl + \OL2\ 

where 



fr/.N rr /, 4ai|a 2 |sin 2 7_t //^l - M V \ cj 2 ± , in v 

u V = u f~ {ai + \a 2 \y ' ^ ) = arctan U^TMJ 7 J' 7± = ^~ • (10) 

As discussed in more detail in this solution describes an exchange of energy between waves 
travelling in the two components u and w. The amplitude of the wave propagating in the u component 
varies between its initial value U, and a smaller non-zero value Umi n = U \ai — |o; 2 | \ / (ai + la^D- There 
is a full exchange in the wave amplitudes from zero to non-zero values when 

(c 2 - l)k 2 = S 2 - 1. (11) 

Then v\ — v\ = v 2 , oli = — a 2 = and the solution ® reduces to 

u = U cos 7_t cos(kx — 7+t), w = — E7 sin 7_t sin(/cx — 7+t). (12) 

This is illustrated in Figure 1, for the case c = 5 = l,k = 1.6. If c = 1, then ([TT]) can be satisfied 
only for 5 = 1, while & is arbitrary. If c ^ 1, then cannot be satisfied if (S 2 — l)/(c 2 — 1) < 0. 
Otherwise, (|11[) is satisfied for a single wavenumber k = ((S 2 — l)/(c 2 — l)) 1 ^ 2 . 




Figure 1: Full exchange in the linear two- wave solution (|12j) . 



For the solution JTJJ), the energy is partitioned according to 

where, for these linear solutions, E c is calculated using only the quadratic term of J5J). Using 7 2 +7+ = 
S 2 + /c 2 , one may show from JSJ and (|T3J) that E = U 2 (S 2 + /c 2 ) /2. From ((15)) we see that there is a 
partial exchange of energy between u and with period T = 7r/7_ . However, if 7_ is small enough, i.e. 
if the frequencies ui and c^ 2 are close enough for some wavenumber fc, then at t = 7r/27_, 37r/27_, • • • 
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the ^-component is almost in equilibrium and almost all of its energy is lost, whilst at t = 0, 7r/7_, • • • 
the ^-component is almost in equilibrium and almost all of its energy is lost. If v 2 ^> then the 
period T of the energy exchange tends to infinity. 

There exists a corresponding solution to ® consisting of a superposition of two left-propagating 
waves. A combination of this left-propagating solution and the right-propagating solution (0 is 
obtained by using the initial conditions 

u(x,0) = Ucoskx, u t (x,0)=0, w(x,0)=0, w t (x,0)=0. (14) 

We then obtain a standing wave involving two pairs of counter-propagating waves: 

u = U(t) cos(7 + t + 6(t)) cos kx, w = ^ ai \ a ^\U g - n ^ cos/cx, (15) 

OLl + \OL2\ 

with U(t), 0(t) and j± as given by (fTT)j). The amplitude of the standing wave in the u component 
varies between U and U m i ni so that once again there is only a partial exchange between the amplitude 
of the components, unless condition (|11[) is satisfied. Then, we have a full exchange, and (|15[) becomes 

u = U cos j-t cos 7 + £ cos/cx, w = U S~ r sin 7_£ sin 7 + £ cos/cx, (16) 

where we note that 7_ < 7+, from (fTU|) . This standing wave solution is illustrated in Figure [2j for the 
case c = S = 1, k = 1.6. 




The energy in solution (|TB|) is E 1 = U 2 (S 2 + /c 2 )/4, although its partition is rather more complicated 
than (|13j). Nevertheless, the general pattern of a partial energy exchange on a timescale 7r/7_ remains. 
For instance, for 7_ < 1 we can say that u loses almost all of its energy near t = 7r/27_, 37r/27_, • • • 
etc. Some examples of this exchange will be given in the numerical simulations of Sectional 

Note that the role of t and x is interchangeable here; one can consider a superposition of waves 
with one and the same frequency, but different wave numbers. 



2.2 Nonlinear solutions 

An exact nonlinear two- wave solution describing energy exchange between the two components u and 
w of the system can be found if c = 1, i.e., if the characteristic velocities of the components 
coincide (2Hj. It is given in terms of Jacobi elliptic functions: 

(w) = 1 + 5 2 ( arcsin (1) + arcsin (j) 2 > ^ 2 = K sn (^ x vi,2(k)t + O , «), (17) 



where o;i = fc, CJ2 = Vl + £ 2 + & 2 , from J7J). When £ 2 = 1 (i.e. z^ 2 = z^ 2 )? the nonlinear solution ([T7|) 
describes a full periodic exchange between the two components of the system (see [28] for details). 
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The form of the solution for small amplitude can be calculated by taking the limit k — > in (|17|). 
Using the well-known approximation (e.g. 

sn(x, k) = sinx — -k 2 (x — sin x cos x) cosx + 0(ft 4 ), 

one can eliminate the secular terms by renormalizing the frequencies (k) in the solution (|17|). and 
thus find an asymptotic expansion with leading order terms: 

t)=Ti^( sin ( te --^ + ^) (;) 

+ sin(/cx-^ 2 (^ + ^ + ^2^ 2 t + O(^ 4 )) (li)) +°(^ 3 )' ( 18 ) 

where ^2 = (1 + £ 2 )/4o;2. Setting n = U/2 and #0 = 7r/2, the leading order terms of (|T%|) reduce to 
the linear solution Q with c = 1. 

In the general case, when c ^ 1, and for a greater number of waves, such exact solutions are not 
generally available. Instead, we study weakly nonlinear solutions by deriving appropriate asymptotic 
models. We do this in the next section, and return to these particular energy exchange solutions later. 



3 Models for weakly nonlinear wave packets 

We consider the weakly nonlinear evolution of pairs of counter-propagating plane waves within system 
To do this we derive equations for the evolution of wave packets, each with the same dominant 
wavenumber k. We suppose that the wave amplitudes are characterised by a small parameter e <C 1, 
and thus introduce the asymptotic multiple-scales expansion 



^3 
w 3 



+ 0(£ 4 ), 



where 



Q e i(kx-uj 2 t) _|_ jj e i(kx-\-iAj 2 t) 



C.C. 



(19) 



The coupled Sine- Gordon equations (0) may then be reduced to coupled equations for the slow spatio- 
temporal evolution of the wave amplitudes A, £?, C and D, which are functions of slow variables. We 
will consider modulations on a long lengthscale x or on a super- long lengthscale X, with corresponding 
slow and super-slow timescales r and T, defined by 



X = ex, X 



: e 2 x, 



et, T 



e 2 t. 



(20) 



We first derive an asymptotic model for modulations on the super-long lengthscale, leading to a system 
of four non-dispersive nonlinear ly coupled envelope equations. We then do the same for modulations 
on the long lengthscale, generally leading to a system of four non-locally coupled NLS equations 
("averaged equations"). We finally present a third model, the coupled NLS equations, which may be 
regarded as a model for modulations on the long lengthscale with both the 0(e 2 ) transport terms and 
the 0(e 3 ) linear dispersive terms retained. 



3.1 The non-dispersive nonlinear ly coupled envelope equations 

To derive this model, one has to suppose that A, B, C and D are functions of X = e 2 x, T = e 2 t (see, 
for example, [1^1 ESI). At leading order, i.e. 0(e) in this expansion, we find a linear wave equation for 
u\ and wi, which is taken to be the four- wave solution ([TU|) . with A, C and D as yet undetermined 
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functions of X, T. At the next order, one obtains solvability conditions in the usual way from the 
equations for u 2 and w 2 , yielding 

A T + v gi A x = m{\A\ 2 + 2\B\ 2 )A + in 2 (\C\ 2 + \D\ 2 )A + ifi 2 B*CD, 
B T ~ v gi B x = -itn{\B\ 2 + 2\A\ 2 )B - ifi 2 (\C\ 2 + \D\ 2 )B - i/i 2 A*CD, 
C T + v 92 C x = m(\C\ 2 + 2\D\ 2 )C + in4(\A\ 2 + \B\ 2 )C + i^D*AB, 
D T ~ v 92 D x = -itis{\D\ 2 + 2\C\ 2 )D - im(\A\ 2 + \B\ 2 )D - i^ A C*AB, 

where the group velocities v 9i are given by 

k 1 + a 2 S 2 c 2 



(21) 



w< l + a 2 S 2 



(22) 



and the nonlinear coefficients are 



S 2 (l -a x f o 2 (l- ai ) 2 (l-a 2 ) 2 



Ml — 1 rr~ TFr\> 



5 2 (l-a 2 ) 4 5 2 {l- ai ) 2 (l-a 2 ) 2 



(23) 



4^ 2 (l + a^ 2 )' ^ 2a; 2 (l + al5 2 ) 



The equations (|2T} are non-dispersive. 

Note that for the case of two-wave interaction, with B = D = for instance, (|2T]l yields the system 

Ar + ^Ax = i{^\A\ 2 + ^ 2 \C\ 2 )A, C T + v g2 C x = i(^\C\ 2 + f i A \A\ 2 )C, 

which is integrable in quadratures ^01 ( see a l so EHD- 

3.2 The non- locally coupled NLS equations 

To derive the second model, we generalise the approach of 16 to the case of four waves. The 
generalization is not straightforward, since the number of characteristic variables now differs from the 
number of independent variables, whereas they are the same in the two- wave case. However, these 
potential difficulties can be readily overcome, and the approach we describe here can generally be 
applied to any number of waves. 

We suppose that A, £?, C and D are functions of the long lengthscale x — ex i m which case they 
must also depend on both slow timescales r = et and T = e 2 t. Once again, at 0(e) (i.e. the leading 
order in this expansion), we find a linear wave equation for u\ and w±, which is taken to be the four- 
wave solution (p"9j) . with A, £?, C and D as yet undetermined functions of r and T. At 0(e 2 ), there 
are solvability conditions A T + v 91 A x = 0, B T — v 91 B x = 0, C T + v 92 C x = 0, and D T — v 92 D x = 0. 
Introducing the variables 

Vi = X~ v 9i r, & = x + v 9i r, (24) 

the solvability conditions give 

A = A( m ,T), B = Bfa,T), C = C( m ,T), D = D(£ 2 ,T), (25) 

i.e. there is a simple functional dependence for each wave amplitude in its own reference frame. If 
we suppose that our solutions have spatial period I in x, which must satisfy I ~ e~ x for modulations 
on the long lengthscale, then A, £>, C and D are periodic in 771, £1, 772 and £2 respectively, with period 
A = el. 

The 0(e 2 ) terms in the expansion may then be written as 

C 2 (m, r, jy***-" 2 *) + D 2 (6, r, T)e < ( fca!+a;2t >l f 1 ) + c.c. (26) 
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At 0(£ 3 ), one obtains solvability conditions in the usual way: 



d T A 2 +A T = -u f lA mm + i^(\A\ 2 + 2\B\ 2 )A + *> 2 (|C| 2 + \D\ 2 )A + i» 2 B*CD, (27) 



z 

|2 , ol /I |2\ 



d T £? 2 + Bt = —t^B^ - i»i(\B\ 2 + 2\A\ 2 )B - i^QCf + \D\ Z )B - i^ 2 A*CD, (28) 
d T C 2 +C T = l -w 2 C mm + i^(\C\ 2 + 2\D\ 2 )C + ifi 4 {\A\ 2 + \B\ 2 )C + ifi 4 D*AB, (29) 
d T D 2 + D T = - l -io'J,D i2i2 - i^{\D\ 2 + 2\C\ 2 )D - i^{\A\ 2 + \B\ 2 )D - ifi 4 C*AB. (30) 

Here the partial derivatives of A 2 , B 2 , C 2 and D 2 with respect to r are evaluated at constant 771, £1, 
772 and £2 respectively, consistent with their functional form ffify. the nonlinear coefficients are defined 
as before in lj2"3jl . and 

„ 1 - t; g 2 + a 2 o 2 (c 2 - ^|) + 4oj 2 a i (v gi - k/u i ){v gi - c 2 k/uj) 
Ui ~ Wi (l+a,?J2) ' (31) 

We proceed further by integrating (|27|) with respect to r at constant 771 from to tq, and dividing by 
To- Letting tq — > 00, and demanding that A 2 remains bounded in this limit, yields 



A T = l -uo f iA riiril +iii 1 \A\ 2 A + iA lim (— [ ° {2 f i 1 \B\ 2 + /i 2 |C| 2 + /i 2 |#| 2 ) dr 

2 ro^oo ^r Jo 



z> 2 lim ( — / B*CDdr ) , 



where we have recalled that A = A(rji,T), and where the integrals are taken at constant 771. Noting 
the functional forms (|25|). we may use ((23)1 to convert the terms involving |£?| 2 , |C| 2 and |D| 2 into 
integrals with respect to £1, 772 and £2 respectively, still at fixed 771, eventually yielding 

A T = ^u; , ;A riim +iii 1 \A\ 2 A + ifi2 lim (— [ ° B*CDdr) 

2 ro^oo \To J J 

1 / fni+To prji+ro rVi+ T o \ 

+ a lim - 2/ii / |S| 2 da+M2/ |C| 2 dr72 + M2/ |£>| 2 d& • (32) 

ro^ooro \ Jr,! Ax J Vl J 

Here ro has been rescaled in each of the final three integrals, and to perform this rescaling we have 
assumed that v gi and v 92 are both O(l), with (v 91 — v g2 ) = 0(1). Finally, we use the periodicity of 
B, C and D to convert the final three integral terms of l|32|) to averages over one period. Performing 
similar operations for (|2 t SH . H'29(l and (|3U|> gives 

A T = 1 -uj'{A ViVi + ifn (\A\ 2 + 2\B~f) A + i^ 2 (jcf + |£>f ) A + i/j 2 \im^ £° B*CD dr^j , 



Bt = - (2\Af + |B| 2 ) 5 - ifi 2 (jq 2 + Jbf^j B - i^ 2 lirn^ (J- J ° A* CD dr 

C T = l -oJ^C mm + iju 3 (|C| 2 + 2\Df) C + ifi 4 (jAp + ji?f) C + ifi 4 lirn^ (J- ^ ° D*AB dr^j , 

D T = - % -u 2 D^ - ifis (2\Cf + |D| 2 ) D - i^ 4 (JXJ + JbJ ) D - i^ 4 Yltc^ (J- £ C* AB dr^j , 

(33) 

where 

JaJ(t) = ±£ \A\ 2 d m , JbJ(t) = \£ |B| 2 d6, 

]df(T) = \C\ 2 d V2 , ]bJ(T) = \£ \D\ 2 d&. 
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These non-locally coupled NLS equations govern the weakly nonlinear evolution of the wave amplitudes 
in the general case, when v gv v 92 and (v 91 — v 92 ) are all O(l). 

If v 91 , v 92 are 0(1) but (v gi — v Q2 ) = 0(e), then the averaging process described above fails, and 
we write 

v 91 =v g -eA, v 92 = v 9 +eA, where v g = V -M±±^± , A = Ml . (34) 

Then, rather than C3, the solvability conditions at 0(e) give A = A(rj, T),B = B(£, T), C = C(rj, T), 
L> = D(£,T), where 

V = X- V g T : £ = X + ^T. 

At 0(e 3 ) we obtain equations similar to ([2T )) -([5U )) . but with additional terms — AA X , A£? x , AC X , 
—AD X on the left-hand side. We then once again perform an averaging operation, yet this time A 
and C are now both functions of rj and T, whilst B and D are both functions of £ and T. Thus, the 
final integral terms can be considerably simplified, and rather than we obtain 



A T - AA^ 



£ T + A5^ 
C T + ACV 



D T - ADj 

with the coefficients defined as before. There is now a local coupling between the A and C modes, 
and between the B and D modes. 

If one permits v 9l and v 92 to be 0(e), then three other scenarios arise in which the averaging 
process fails. If v 91 = 0(e) but v 92 = 0(1), then using the obvious extension of the above analysis 
one obtains a system in which A and B are locally coupled, whilst C and D satisfy non-locally 
coupled NLS equations. If v 92 = 0(e) but v 91 = 0(1), then C and D are locally coupled, whilst A 
and B satisfy non-locally coupled NLS equations. If both v 91 and v 92 are both 0(e), then something 
slightly different occurs. Rather than ([25)1 . the solvability conditions at 0(e) give A = A(x,T) : B = 
B(x,T),C = C(x,T),D = D(x,T). At 0(e 3 ), the averaging operation immediately yields a system 
of coupled NLS equation for A, £?, C and D, with no non-local coupling, which is equivalent to system 
([55)) introduced below. However, none of these scenarios are further considered here, since v 91 and 
v 92 are typically 0(1) for the Sine-Gordon system J5J). 

3.3 The coupled NLS equations 

An alternative approach is to seek a system with linear dispersive terms, unlike ([21)1 . only local 
nonlinear coupling, unlike ((55)) and ([55)). and with no a priori restrictions on the size of the group 
velocities, unlike ([55)). ([55)) and the other systems mentioned at the end of Section 15^1 Such a model 
is obtained, for instance, by allowing the wave amplitudes A, £?, C and D to depend on the long 
lengthscale x — £X an d the slow timescale r = et, whilst retaining all terms up to order 0(e 3 ). The 
result is 



A T + 


v gi A x - 


1 II A 

2 £u> i A xx 


= i f i 1 e(\A\ 2 + 2\B\ 2 )A + 


i^e{\C\ 2 + 


\D\ 2 )A + in 2 eB*CD, 




v gi B x + 


-eu"B xx 


= -ifi 1 e(2\A\ 2 + \B\ 2 )B 


-i(i2e(\C\ 2 


+ \D\ 2 )B -ifi 2 sA*CD, 


C T + 


- V 92 C X - 


-eu>' 2 ' C xx 


= i f i 3 e(\C\ 2 + 2\D\ 2 )C + 


i/j,4e(\A\ 2 + 


\B\ 2 )C + ifi 4 sD*AB, 




V 92 D X + 


2 £UJ 2 D xx 


= -i» 3 e(2\C\ 2 + \D\ 2 )D 


-ili 4 e(\A\ 2 


+ \B\ 2 )D - i/i A eC*AB. 



(36) 



= % -J[A m + ifi, (\A\ 2 + 2|B| 2 ) A + ifi 2 (|Cf + |D| 2 ) A + i^i 2 B*D C, 

= - V% - m (2\Aj+ \B\ 2 ) B - i M2 (\C? + \D\ 2 ) B-ifi 2 A*C D, 

. 2 V 1 V 1 (35) 

= l^'Cw + ^3 (|Cf + 2\Df) C + m (\A\ 2 + \B\ 2 ) C + D*B A, 

= " ^ 3 ( 2 K^+ l^l 2 ) D ~ *M4 (W + \B\ 2 ) D-ifi 4 C^A B, 
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In this form the linear dispersive and nonlinear terms are of the same relative order, but 0(e) smaller 
than the transport terms. If instead we had used the variables X = e 2 x, T = £ 2 t, then we would 
have had a balance between the transport terms and the nonlinear terms, with the linear dispersive 
terms being relatively 0{e 2 ). Only in certain special cases may one simply rescale x and t so that all 
terms of (|36|) are simultaneously of the same order. For instance, if only one wave is present, then in a 
reference frame moving at the group velocity of the wave one may use T = e 2 t and x — ex - Similarly, 
if B = D = 0, and v 91 = v 92 , one may perform the same operation to obtain coupled NLS equations 
for A and C. 

Nevertheless, is an appealing model on purely physical grounds. Higher order terms, both 
linear and nonlinear, could, in principle, be dynamically important. Indeed, since the effects of 
transport, dispersion, and leading order nonlinearity are captured, models of this type are widely used 
in studies of weakly nonlinear wave phenomena, as mentioned already in the Introduction. There is an 
analogy here with the well-known Boussinesq equations often used to model the two-way propagation 
of water waves. Indeed, any of the asymptotic models and (|35|) already introduced, or 

those mentioned at the end of Section ET^l can be derived from by an asymptotic multiple-scales 
reduction on the appropriate length and time-scales. 

Note that one can deduce (|36|) without recourse to a new asymptotic expansion. With the usual 
approach (e.g. JSHEEI)? the terms on the left-hand side of (|36|) result from a Taylor series expansion 
of the dispersion relation about the dominant wavenumber of that packet. Only the first two terms 
are retained. The terms on the right-hand side are the relevant nonlinear terms, as already obtained 
in (|21|). for instance. 



4 Modulational instability of two-wave solutions 

It is well known that dispersive weakly nonlinear systems of the type being studied can support 
modulational instabilities. Some simple examples of this for two- wave solutions were given in our 
previous study [TO] . These results are briefly summarised in Section 14.21 Then, in Sections 14.31 
and 14.41 the analysis is extended to the other asymptotic models derived in Section 03 



4.1 Plane waves 

We consider the modulational instability of spatially uniform coupled plane wave solutions of the form 

A = A e iQaT , C = C e iQcT , B = D = 0. (37) 

This is a solution of the non-dispersive nonlinearly coupled envelope model, the non-locally 
coupled NLS models and and the coupled NLS model provided 

Ox = /il|A | 2 + /i 2 |C | 2 , = /i3|C | 2 + /i 4 |A | 2 . 

Substituting (|37|) into (|T§|) we obtain the leading order solution in the form of two counter-propagating 
waves: 

Q = £Ao S kx -^-^«)t) + ^^(kx-^-e^t) + c c + 0{£ 2y (38) 

This solution describes a weakly nonlinear interaction between the two components. For e ^ 0, it 
generalises the linear energy exchange solutions of Section bv nonlinearity-induced corrections of 
0(s 2 ) to the frequencies uoi : 2(k) of the linear waves, given by Q. For e —> 0, it reduces to the two- wave 
linear solution ® if 

£A ° = 2(aT+ U \a 2 \y £Co= 2{aT+\a 2 \Y (39) 

In the case c = 1, matches the weakly nonlinear expansion (|TH|) of the exact two- wave energy- 
exchange solution (|TTj) . Since a\ = 1 and (12 = —S~ 2 at c = 1, from ((23) Mi = = /M = 0, and 
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M3 = (1 + S 2 ) 3 /(4u 2 S 4 ), so that taking 

we find ft a = 0, s 2 Q c = /i3|eCo| 2 = ft 2 (l + ^ 2 )/4cj2, in agreement with ([T8J). 
4.2 Previous results 

In [TUI we considered the stability of the coupled plane wave solutions ((57)) in the context of two of 
the asymptotic models of Section 03 Firstly, we showed that the non-dispersive nonlinearly coupled 
envelope equations (|2U with B = D = do not support any modulational instabilities. Secondly, we 
studied instabilities with modulations on the long lengthscale x = ex -> w ^h v 91 — v 92 non-zero and 
O(l), so that we used ((33)) with B = D = 0. However, the dynamics of that system are relatively 
simple, since as in [T7| the equations may be transformed to two uncoupled NLS equations using the 
nonlocal transformation 



A( m ,T) = A(r h ,T)exp jf \C\ 2 (T')dT^j , 
C{ m , T) = C{ m , T) exp [i^ £]Af (T')dT^ . 



(40) 



Alternatively, stability results can be established directly by looking for solutions of of the form 
i4(jh, T) = (l + aie i (^ 1 -"" T ) + a* 2 e-< km -n:T)^ , 

C( m , T) = C (T) (l + de^ 2 ""^) + c* 2 e-< kT >2-n* c T)^ ^ 



(41) 



representing sinusoidal perturbations to the coupled- wave solution (|57j). Neglecting terms quadratic 
in disturbance amplitude, (|5Sj) leads to two uncoupled linear systems, as expected from (^UJ), with 
corresponding dispersion relations 

^ = - Ml |A | 2 ^ 2 + ^, ^ = -^|C |M'« 2 + ^- (42) 

Using ((201) an d lEH we can rewrite these in terms of an unsealed disturbance wavenumber /c and 
frequency a), via 

sk = k, e 2 Cl a — uj a — kv gv s 2 £l c = lu c — kv Q2 , (43) 

giving 

(u) a -kv g ^ =-/ii|^ |V/P + ^^, (tD c -fcv P2 ) = -^skC^ol 2 ^^ 2 + — , (44) 

where for formal validity we need Ao, Co = 0(1), k = O(e), and — v 92 ) = 0(1). 

From ((25)) . we see that /ii > and /i3 > 0, so there will be instability only if ou" > or ou 2 ' > 0. 
Using JSU), we can show that o;^ is always positive. Thus for sufficiently small k there will be a root 
for ou c with positive imaginary part, corresponding to an unstable mode. Using 0, we can show that 
ou" < for small fc, whilst ou" > for large /c, tending to zero at infinity. Thus ou a is real for small 
fc, and hence does not correspond to an unstable mode, whilst for larger /c, where ou" > 0, it will 
correspond to an unstable mode for sufficiently small k. 
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4.3 New results: the coupled NLS equations 

We now consider instability in the context of the two other asymptotic models of Section 01 The first 
model is for modulations on the long lengthscale \ — ex i but with (v 91 — v 92 ) = 0(e). Then we must 
use equations with B = D = 0, giving a system that is locally coupled in A and C. The second 
model is the coupled NLS equations with B = D = 0. However, with the appropriate rescaling, 
it is easy to see that these two systems of equations are equivalent. Therefore, even though the two 
models were derived in different ways, they may be analysed together. 
We look for disturbances to the coupled- wave solution (|37|) in the form 

A = A (T) (l + ai e^ kx -^ + aJe-^-^c*)) ; 

C = Co(T) (l + c^Ckx-Zact) + c * e -i(kx-u;: c t)^ ^ 

where we find it convenient to revert to the original unsealed variables x, t and work in terms of the 
unsealed disturbance frequency u and wavenumber k. Substituting into or ((36)) we obtain: 

u ac ~ kv 91 - ^u"k 2 + /ni\eA \ 2 ^j ai + /ni\€Ao\ 2 a 2 + /i 2 |£C | 2 ci + /i 2 |£C | 2 c 2 = 0, 

1 

fi4,\sA \ 2 ai + /i4Mo| 2 «2 + (^) ac - kv 92 - ^u 2 k 2 + ^ 3 |eC | 2 J ci + /i3|eCb| 2 C2 = 0, 

-/i4\sAo\ 2 ai - /J J A\sA \ 2 a 2 - /i 3 |£Co| 2 ci + ^a) ac - fcv P2 + ^ 2 '^ 2 ~ M3kCo| 2 ^ c 2 = 0. 

This four-by-four linear homogeneous system of equations must have a vanishing determinant, from 
which we find that 

Qa(v, k)Q c (v, k) = R, where R = \eA \ 2 \eC \ 2 , (45) 



-/ii|£A | 2 ai + ( cu ac - kv 91 + -u"k 2 - m\sAo\ 2 ) a 2 - HeCopci - /i 2 kC | 2 c 2 = 0, 



Q a (fi, fc) = (t> - ^) 2 - + /W/Mo| 2 , Qc(«, *) = («- ^ 2 ) 2 - + W^kCol 2 - (46) 

Here v = a) ac / fc is the phase speed of the disturbance, whilst R is proportional to the coupling constants 
in (|35|) and (|36|) . Note that Q a = and Q c = are the dispersion relations for the uncoupled NLS 
equations for A and C respectively, cf. (|^|) . Equation is a fourth-order dispersion relation with 
real coefficients for v. Complex solutions for v imply modulational instability. 

4.4 Instability analysis: the coupled NLS equations 

Working within the context of (|35[). for formal validity of (|45|) we need Ao, Co = 0(1), k — 0(e) ^ and 
(v gi — v 92 ) = O(e). However, within the context of there are no specific scale restrictions, and 
we may investigate different scaling regimes. 

We first probe the effects of the dispersive terms for modulations on the super-long lengthscale, by 
taking k = e 2 K, with K = 0(1) and e« 1. Writing v = i)\ + e 2 v 2 + • • • , one may derive expressions 
for the four roots as 

(fi „ J2 n , tr\ a |2 , 4 A, 7/2^2 ■ /i2/i4^i^ 2 |A Q | 2 |C | 2 \ 6 

/V, „ ^ ^2,, ,//|^|2 , „4 Z" 1 , ,//2fr2 , /i2M4^i^ 2 |^o| 2 |Co| 2 ^ n , 6 v 

V 4 V^l ^2/ / 
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Since /as > and uo' 2 ' > 0, at least one pair of roots for v is complex, and hence there is always an 
unstable mode. Since fi\ > 0, there will be an additional unstable mode if uj'{ > 0. Thus, our analysis 
suggests that the dispersive terms are destabilising for disturbances on the super-long lengthscale, 
leading to an instability with a growth rate of 0(e 3 ). This is consistent with the nondispersive 
analysis of mentioned in Section FOl which showed there will be no instabilities on the super- 
long lengthscale with growth rates of 0(e 2 ). 

We can recover results for modulations on the long lengthscale by taking k = ek, with k = 0(1) 

and £ <C 1. If (v 91 — v 92 ) = 0(1), then writing v = i)\ + ev 2 H , one may derive expressions for the 

four roots as 



„2, 2 \ , £V2/"4wi'w 2 'l^0| 2 |C | 2 +0 ( £ 6) 



In r Y - - 2 / „ , tn r |2 !, //2 £ .2\ , £ 4 A*2A*4^l ^2 1^0 1 2 1 ^0 P , n( c 6\ 

\ v - v a2) - ~ £ M3^ 2 l G o| - t w 2 K H 7 —^2 VU(e). 



( V 91 - V 9 2 ) 2 
( V 91 - V 9lf 



(47) 



Since {13 > and uj 2 > 0, there will be an instability with a growth rate of 0(e 2 ) if k 2 < 4fis\Co\ 2 /uj 2 . 
Since /ii > 0, there will be an additional unstable mode if ou" > and k 2 < 4fii\Ao\ 2 /ou" . This is 
consistent with (J12|) and (@3|), derived formally from the non- locally coupled NLS equations. If (v 91 — 
v 92 ) = O(e), then no further approximation of ((45)) is possible, since it is already the asymptotically 
consistent model for this regime. 

More generally, we now consider the properties of the roots v(k) of (|35)> without any scale restric- 
tions on or e. First, if k ^> 1, then the roots of (|l5|) are approximately £ = ±uu"k/2 and # = ±0;^ fc/2, 
so that there is stability. We then consider how the roots vary as k is reduced from infinity, with all 
other parameters fixed. Note that these other parameters consist of the plane wave parameters /c, 
\sAq\ and |sCo|, and the system parameters c and 5. We consider the graph of Q = Q a (v, k)Q c (v, k) 
as a fourth order polynomial in as k is varied. Stability occurs whenever this graph intersects the 
line Q = R four times, and otherwise there is instability. Note that R is independent of k and so this 
line is immovable as k is varied. Noting that /i2, /M and ou 2 are all positive, R takes the same sign as 
ou". There are thus two possible cases to consider: 

1. Suppose first that ou" < 0, so that R < 0. For sufficiently large /c, both uncoupled modes are 
stable and so Q a Q c has four real zeros, = 1,2,3,4 where v\ > £2, £3 > £4 are the zeros 
of Qa,Qc respectively. Since fi\ > and ou" < here, are real for all k : from ([ifi]). But 
as fc is reduced from infinity, there is a critical value of k = fc C2 = 2 |sOo 1(^3/^2 ) x ^ 2 , at which 
v 3 = v 4 = v 92 , and below which become complex, so that the uncoupled O-mode is then 
unstable. The issue then is whether the coupling in suppresses or enhances this instability. 
Consider then the situation when k = k C2 , so that Q c = (v — v 92 ) 2 and Q a = (v — v 91 ) 2 — f3 2 , 
where 

(i 2 = KK) (m3<|£C | 2 - m<4'Mo| 2 ) > 0. (48) 
There are two sub-cases to consider: 

(a) If \v Q2 — v 91 \ > (3, then the double zero v Q2 of Q c (v,k C2 ) does not lie in (v 2 ,vi). This 
scenario is depicted in Figure [Hh- Since R < 0, Q a (v, k C2 )Q c (v, k C2 ) = R must have at least 
two complex roots, so that the coupled system is unstable at k = k C21 and therefore for 
some k > k C2 . In this sense we can say that coupling is destabilizing. Indeed, in this case 
the system remains unstable for all k < fc c , where k c > k C2 is the critical value at which 
the system first becomes unstable, and is the value of k at which the line Q = R is tangent 
to the graph of Q = Q a Q c - 

(b) If \v 92 — v 91 \ < /3, then the double zero v 92 of Q c (v,k C2 ) lies in (v2,v\). This scenario is 
depicted in Figure Ob- Stability at k = k C2 is then determined by 



i? c = max(Q a (v ± , k C2 )Q c (v±,k C2 )), v± = v 92 H 



Modulational instability of two pairs of counter- propagating waves in two- component media 15 



(a) R<0 (b) R<0 (c) R>0 (d) R>0 




v2 v- v1 v+ vg2 v2 v- vg2 v+ v1 v2 v- vg2 v+ v1 v2 v- v1 v+ vg2 



Figure 3: Schematic representations of the possible stability scenarios at k = k C2 . Plotted in black in each 
case is the coupled dispersion function Q = Qa(v,k C2 )Qc(v,k C2 ) as a function of v, as defined in @3J- The 
horizontal lines denote values of R with differing stability properties. 

where v± are two of the turning points of Q = Q a {v, k C2 )Q c (v, fc C2 ), the third being at 
v QT When R < R C: Q a (v/k C2 )Q c (v/k C2 ) = R has at least two complex roots, and we 
can say that coupling is destabilizing, whereas if R c < R < there are four real roots of 
Qa(v,k C2 )Q c (v,k C2 ) = R, and we can say that coupling is stabilizing. But here, in both 
scenarios, the stability situation may change again as k is further reduced below k C2 . 

2. Suppose next that uj" > 0, so that R > 0. For sufficiently large fc, both uncoupled modes are 
stable, whilst for sufficiently small k, both uncoupled modes are unstable. Suppose that as k is 
reduced from infinity, one of the uncoupled modes becomes unstable for some value of fc, while 
the other, for this value of k is still stable. To be explicit, let us suppose that it is the C-mode 
which becomes unstable, so that k = k C2 as above. The other case can be treated in the same 
way. Now at k = k C2 our hypothesis that the A- mode is stable again implies that f3 2 > 0, and 
once again there are two sub-cases to consider: 

(a) If \v 92 — v 91 \ < /?, then the double zero v 92 of Q c (v/k C2 ) lies in [^^l]- This scenario is 
depicted in Figure^. Since R > and all turning points of Q = Q a (v,k C2 )Q c (v, k C2 ) lie 
in Q < 0, the coupled system is unstable at k = k C2 , and therefore for some k > k C2 , and 
in this sense we can say that coupling is destabilizing. Indeed, in this case the system 
remains unstable for all k < fc c , where k c > k C2 is the critical value at which the system 
first becomes unstable, and is the value of k at which the line Q = R is tangent to the 
graph of Q = Q a Q c . 

(b) If \v 92 — v 91 \ > (3, then the double zero v 92 of Q c (v/k C2 ) does not lie in [#2,£i]- This 
scenario is depicted in Figure EJi. Stability at k = k C2 is then once again determined by 
R C1 defined as before, which is now the value of the single turning point with Q a Q c > 0. 
When R > R Cl Q a (v, k C2 )Q c (v, k C2 ) = R has two complex roots so coupling is destabilizing, 
whereas if < R < R c there are four real roots of Q a (v, k C2 )Q c (v, k C2 ) = R and coupling 
is stabilizing. But here, in both scenarios, the stability situation may change again as k is 
further reduced below k Co . 



4.5 Summary 

We have considered how modulational instabilities of spatially uniform coupled two-wave solutions 
are predicted by the four weakly nonlinear models of Section |2l As noted in Section FOl the non- 
dispersive nonlinearly coupled envelope equations ([2T]) give an unambiguous prediction of stability for 
disturbances with k of 0(e 2 ) on a timescale ~ e~ 2 . Further, the non- locally coupled NLS equations 
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([55)1 give an unambiguous prediction of instability for disturbances with k of 0(e) on a timescale ~ £ -2 , 
whenever v 91 — v 92 is 0(1), and either uj'{ > or uJ^ > (the latter is always true here). As shown in 
Sections l4.3l and rOl the coupled NLS equations (|35|) and (|36j) give a dispersion relation which is harder 
to analyse, but we showed that it reproduces the results for disturbances with k of 0(e) or 0(e 2 ) in their 
respective regimes of validity. However, in the context of (|36|). the coupled NLS dispersion relation 
also gives predictions which are strictly outside its asymptotic range of validity. Such predictions 
need to be tested independently, for instance, by numerical simulations of the full unapproximated 
Sine-Gordon equations. We perform such tests in Section^ by numerically determining the behaviour 
of system (|5J. 



5 Modulational instability of four-wave solutions 

We now discuss the modulational instability of four- wave solutions. We use the same ideas as in 
Section^ except that we restrict our analysis to the non-locally coupled NLS equations, derived for 
modulations on the long lengthscale \ = ex - 



5.1 Plane waves 

We start by examining the form of spatially uniform coupled plane waves. For such solutions, the 
non-locally coupled NLS models ([55)) and ([55)) reduce to the same system. Looking for solutions of 
the form 

A = A e iQaT , B = Boe~ iQbT , C = C e iQcT , D = Doe~ iQdT , (49) 
we see that a balance is only possible if 

n a -n b = n c - n d . (50) 

Substitution of into (|55J) or ((53J) gives 

V a = /ii(|A | 2 + 2|5 | 2 ) + /i 2 (|C | 2 + |A)| 2 ) + ^/\M\ 
Q b = /ii(2|A | 2 + \B \ 2 ) + /i 2 (|C | 2 + \D \ 2 ) + M20*/I^o| 2 , 
^ c = /i 3 (|C | 2 + 2|A)| 2 ) + /i 4 (|A | 2 + \B \ 2 ) + ^/\C \ 2 , 

n d = /i 3 (2|c | 2 + \d \ 2 ) + M\M 2 + l^ol 2 ) + ^H\Do\\ 

where 

<l> = AoBoC^D*. (51) 
Thus, for solutions of the form (|49|) . the condition ([KU)) for balance becomes 

d^P _ M ( w + e&fi . (|CoP _ |D „ n ( w + ^) . (52) 

If (O is satisfied, then we may substitute (03 into ([T9)) to obtain a solution in the form of two pairs 
of counter-propagating waves, 

+ £ ^Coe^^-^ 2 -^^)^ + jDoe ^(^+(- 2 -- 2 ^)^)) + c.c. + 0(e 2 ). (53) 

This solution generalises the linear energy exchange solutions of Section bv nonlinearly-induced 
corrections of 0(s 2 ) to the frequencies ooi^ik) of the linear waves, given by J7J). For e —> 0, it reduces 
to the four- wave linear solution ([13)) if 

eA = = ' — jr, eC = eD = - ( — — rr, (54) 

4(ai + \a 2 \) 4(ai + |a 2 |) 

which clearly satisfies the balance condition ([52)) . 
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5.2 The non-locally coupled NLS equations: (y 91 — v 92 ) = 0(1) 

We first study the instability of the plane- wave solution in the case with v 91 and v 92 both O(l), 
and (v 91 — v 92 ) non-zero and 0(1). Accordingly we use equations (|33|) . and we look for solutions 
in the form (|41j). with corresponding forms for B and D. Neglecting terms quadratic in disturbance 
amplitude, this yields four linear systems of the form 

(n a - X -J{k 2 + /ii|A | 2 - ai + /ii|A | 2 a 2 = 0, 

-fiMol 2 ^ + (fi a + ^'£ 2 - /ixl^ol 2 + a 2 = 0. 

Note that although this system might appear only to involve |A)|, in fact the terms involving give 
a coupling to the other basic state waves via Bo, Co and Do, unlike the two- wave case. Each system 
has an associated dispersion relation of the form 

+ w (0 ~ r) ha " G w "* 2 + w " Ml|A ° |2 ) G w "* 2 + w ~ " l|A ° |2 ) + m2|a ° |4 = °" 

When Aq, Bo, Co and Do all have the same phase, as is the case for ((54)) for instance, then (j) = 
\AoBoCoDo\, and the dispersion relation reduces to 

6 2 f 1 ,,"^ „ I J |2 , Mgo^Dgol V 9U ,4 /rrN 

a = \2 1 ~^il^o| + J - V>i\Ao\ • (55) 

Reverting to dimensional variables using ((43)) gives 

(d, -fc Vgl ) 2 = Q w ;'P-(l-r ) Ml | £ Ao| 2 ) -/xfleA)! 4 , ra = ^ 

with corresponding forms for cD c and a)^. Although equations ([55)) and ([56)) have a similar form 
to their two- wave counterparts (|42|) and ((44j) . there are important differences. For instance, consider 
the instability of the A- mode, as given by 

When a;" > 0, (@3J) shows that for the two- wave solution a maximum growth rate of s* = fii\eAo\ 2 
is achieved by the A- mode when fc 2 = fc 2 = 2/ii|£Ao| 2 /|^i | > 0. However, for the four- wave solution, 
instability of the A- mode is instead determined by the parameter r a , and only if r a = are the 
two- wave results recovered. If < r a < 1, then ([56|) is minimised at fc 2 = fc 2 (l — r a ), so that the 
disturbance wavenumber of the most unstable A-mode is reduced, although the growth rate remains 
s*. If r a > 1, ([56)) is instead minimised at k = 0, with value s2(r 2 — 2r a ) > — s 2 . Therefore, if 
1 < r a < 2 the most unstable A- mode occurs at fc = with a reduced growth rate s < s*, whilst if 
r a > 2 the A- mode is stable. Thus, a two- wave solution can yield an instability in the A-mode whilst 
a corresponding four- wave solution does not. 

When uj'I < 0, shows that the two- wave solution does not lead to an instability in the A- mode. 
However, when r a > 1, ((56)) is minimised at fc 2 = (r a — l)fc 2 , leading to an instability in the A- mode 
with growth rate s*. If < r a < 1, ([56|) is instead minimised at fc = 0, with value ^(r 2 — 2r a ) > — s 2 , 
leading to an instability in the A- mode with a reduced growth rate s < s*. Thus, a four- wave solution 
can yield an instability in the A- mode whilst a corresponding two- wave solution does not. 

There are corresponding conclusions for the B mode, and for the C and D-modes (for which 
UJ2 > 0), and the overall stability of (|H?)) will be determined by the most unstable of these four modes. 
Although the possible stabilising and destabilising effects outlined above are of interest, we do not 
consider the details any further here. 



BqCqDo 



(56) 
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5.3 The non-locally coupled NLS equations: (v 91 — v 92 ) = 0(e) 

We now turn to the case with v 91 and v 92 both O(l), but (v 91 — v 92 ) = 0(e). Accordingly we use 
equations (|35|) . and we look for solutions in the form 

A(rj,T) = A (T)(1 + cue^-VacT) + a * e -*(^-^)) 5 

C(rj,T) = C (T)(1 + cie^^-^ acT) + c£e- <( ^- A " T >), 

with corresponding forms for B and D. Neglecting terms quadratic in disturbance amplitude leaves 
a coupled system in A and C, and a coupled system in B and D. The latter system reduces to the 
former with the equivalence A^B^C^D^^ rj-,&bd ^ -^ fl c 5 so it is sufficient to just consider 
the behaviour of the coupled A, C system, say, which reduces to solving 

1 , Ma.2 i , . \ a 1 2 M2</>* \ , , , \ a |2 _ , (\n |2 , ^* A . , 1^ |2, 



^ ac + - + /ii|A | 2 - j^|2 J a i + Mil^oi a2 + ^2 (jc r + ) ci + /i 2 |Corc 2 = o, 

-/ii|A | 2 ai + ^ ac + + £ 2 - /ii|A | 2 + jj^j a 2 - /i 2 |C | 2 ci - fi 2 (jC | 2 + ) c 2 = 0, 



| 2 + ^ja 1+/ M|A | 2 a 2 ^^ ■ -^-' 2 ' 



-M4|^orai - ^4 (JA)T + J «2 - MsiCorci + [n ac - kA + - /i 3 |c r + J c 2 = o. 

Thus, as in Section FQ1 the dispersion relation may be obtained by setting a 4 x 4 determinant to be 
zero. However, in contrast to Section ^31 the presence of the terms involving <\> makes it hard to derive 
analytically a quartic equation for tl ac . Thus, we find the roots numerically for a given parameter set, 
by calling an eigenvalue routine for the above system, giving Cl ac directly. 



6 Numerical simulations 

We perform numerical simulations of the original system of coupled Sine- Gordon equations to 
examine the accuracy and validity of the theoretical predictions of instability derived in Sections 0| 
and [5] In our previous study [TUj we verified that the modulational instability of two- wave solutions 
is well described by the non-locally coupled NLS equations, over a wide parameter range. Here, we 
make a more detailed comparison between the predictions of this theory, those of the coupled NLS 
equations, and the numerical results, for the instability of both two- wave and four- wave solutions. 
Although our emphasis will be on verifying the growth rate of modulational instabilities, we will also 
comment on the structure of the solutions when any instabilities reach maximum amplitude. 

6.1 The numerical method 

Since the plane wave solutions and the modulational instabilities we analyse are spatially periodic, it 
is natural to use a Fourier basis for the x-dependence of u and w. Then, denoting a quantity q(x,t) 
as 2j>o Qj(t) e ^ x + c - c -? where j is one of the numerically resolved wavenumbers, we write (|3J) as the 
first order system 
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(57) 



where the terms Nj(t) originate from N(x, t) = (u — w) — sm(u — w). We solve this using a split-step 
scheme, a method which has frequently been applied to such nonlinear wave equations. Neglecting the 
nonlinear terms in ((57)). i.e. those involving TV, we obtain a linear system that can be updated exactly, 
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since given u, u t: w and w t at time £, and uii^(j), we can write down explicitly the solution of the 
linear system at time t + At. Similarly, neglecting the linear terms in JSZJ), i.e. those not involving TV, 
we obtain a system in which u and w are constant, so that Ut and Wt can be trivially updated by using 
the value of N at time t. This latter property relies on the second derivatives with respect to t in the 
governing equations. By combining these two steps in a sequence of appropriately chosen fractional 
steps, one may obtain a method of any order (as explained, for instance, in section 2 of 39 J. For the 
simulations to be presented, we use a fourth order scheme with a time-step in the range 0.025 to 0.1, 
with 96 to 768 points in the ^-direction, and the term TV evaluated pseudospectrally. The nonlinear 
scheme then conserves energy to a high degree - typically to within 10 -6 %. If the nonlinear terms 
are neglected, then the scheme exactly resolves linear dynamics, within roundoff error. 

For initial conditions, we take either (|HJ) or (|T^ which lead to the two- wave or four- wave exchange 
solutions (|39|) and (|5^j) respectively, with wavenumber k and amplitude parameter U. The domain 
width / is then restricted by the need to choose / = 2irn/k, where n is a positive integer. To trigger 
any instabilities, we add background noise to u and w, at a level at least 1000 times smaller than the 
amplitude parameter U. In unstable cases, we identify a modulational instability of wavenumber k 
by calculating the energy in pairs of sideband modes with wavenumber s k±k. To do this we consider 
the quantity 

E side \ m (t) — E k+k + E k _y k = 27rra/Z, with m a positive integer, (58) 

corresponding to the energy in the ra-th pair of numerically resolved sideband modes. Here Ej is 
the energy, calculated according to (0), when u and w are replaced by their Fourier components 
with wavenumber j. Initially, E side \ m (t) typically grows exponentially for one or more m, permitting 
growth rates to be calculated for the corresponding wavenumbers k = 2i\mjl. 

To check the relevance of these growth rates, we also distinguish a disturbance energy. Although 
the linear solution projects solely onto wavenumber /c, the unperturbed nonlinear plane wave solution 
under consideration projects onto wavenumbers (2n + l)/c, for all non-negative integers n. Therefore 
we define an energy associated with the plane wave and its nonlinear harmonics 

Eplane = ^ E( 2n +l)k, (59) 
n>0 

and an energy associated with the disturbance field 

E= Ej. (60) 

j/(2n+l)fc 

An overall disturbance growth rate can then be calculated by analysing E. Typically, the growth rate 
of one pair of sideband modes approximately equals the growth rate of E, so that we can identify with 
some certainty which sideband pair yields the most relevant modulational instability. 

This method has been tested by calculating the growth rate for the modulational instability of a 
one- wave solution at c = 1 and e <C 1, when stability is determined by that of a corresponding uncou- 
pled Sine-Gordon equation. Calculated growth rates typically agree with the well-known theoretical 
prediction (e.g. jSZI) to at least four significant figures. 

6.2 Instability at small amplitude 

We first consider the instability of small amplitude plane- wave solutions, for which the individual 
wave amplitudes eAo, eBq, eCo and sDq are all less than 0.1. We reduce the parameter regime 
to manageable proportions by only considering waves that satisfy the full-exchange condition ([TT|) . 
Then, according to (|39|) and (|54|h and recalling from Section [21 that a± = — ol^ = for full exchange 
parameters, we have eAq = eCo = U/4 for two- wave solutions, and sAq = sBq = eCo = eDq = U/S 
for four-wave solutions. Here U is the maximum amplitude of u in the linear solution, although the 
maximum amplitude of w is U/5. 
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The results are summarised in Table [3 Numerically determined growth rates s n are given for 
modulational instability at just a single wavenumber fc n , corresponding to a relatively strong instability 
in each case. For the two- wave runs U has been chosen so that the maximum value of u or w in the 
linear solution is 0.2, so that the values of sAo and eCo are 0.025 or 0.05. For the four-wave runs, U 
is twice that of the corresponding two-wave simulation, so that the maximum value of u or w in the 
linear solutions is now 0.4. However, the individual wave amplitudes are the same in both cases, and 
in this sense the runs are closely related. 
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Table 1: Numerically determined growth rates s n for small amplitude solutions. 



We first discuss the two- wave results. There are two theoretical predictions for the growth rate. The 
first comes from the non-locally coupled NLS equations, valid when k = 0(e) and (v 91 — v 92 ) = 0(1). 
This theory predicts the existence of two possible modes of instability, one evident as growth in wave 
amplitude A and the other as growth in C, and can be diagnosed from (|44|). Since uo" < for all the 
parameters considered here, Co a is real and hence A is always stable. Thus, we need only calculate the 
growth rate of C, which we denote by s c = lm(cj c ). The second prediction comes from the coupled 
system (|45|) . valid when k = O(e) and (v 91 — v 92 ) — 0(e) , or equivalently from the coupled NLS 
equations. This predicts growth of a coupled mode in A and C, and we denote this growth rate by 
s ac Im(o> ac ). 

As shown in Table d both s c and s ac are in good agreement with s n for all the parameters 
considered. Indeed, it is perhaps surprising that s c is a good approximation for the final one or two 
cases, when v 91 — v 92 is apparently small. However, one can also see that s c and s ac are almost equal 
for these cases. This seems to be because for a wide range of full exchange parameters, such as those 
in Tabled at least one of /i2, an d ^2 * s small. Thus, for small amplitude conditions, the 

coupling term R in is not important, so that even when v 91 — v 92 is small the growth rate is 
approximately determined by the uncoupled instability in C. 

Turning to the four- wave results, we once again have two predictions for the growth rate. The 
first is the uncoupled growth rate for the C and D modes, which have the same growth rate s c . Even 
though s a = 55 can be positive, s c > s a for all the parameters in Table ^ The second prediction, 
denoted by s aCl is no longer derived from the coupled NLS equations, but rather from the non-locally 
coupled NLS equations, valid when k = O(e), with (v 91 — v 92 ) = O(e), and v 91 ,v 92 both 0(1), as 
described in Section lOl It now describes a coupled instability in A and C, or an equivalent instability 
in B and D, with the same growth rate. 

As shown in Tabled there is now generally a small discrepancy between the numerically determined 
growth rate and the best theoretical prediction. We speculate that this is an effect of the wave 
amplitude - recall that in the linear solution the maximum amplitude reached by u or w is 0.4, twice 
that of the two-wave simulations. This point will be clarified in the next section. 
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6.3 Instability at moderate amplitude 

We now consider how s n and the functions s c , s ac change as U is increased, for three different parameter 
sets. We first discuss a full-exchange case with c = £ = 1.5, fc = 1, for which v 91 = 1.23, v 92 = 0.75, 
so that v gi — v g<2 = 0.48. The nonlinear interaction coefficients are fi\ = 0.003, H2 = 0.13, [13 = 1.00, 
/i4 = 0.08, and uj'[ = —0.31, uo 2 = 0.73. The results are shown in Figure 01 For the two-wave case, 
given in Figure St, the theoretical predictions s c and s ac are almost indistinguishable. For small U 
they closely follow s n , whilst as U is increased the differences between the predictions and s n grow. 
The four- wave results, given in Figure^, show the same type of behaviour. However, note how the 
four- wave predictions are different from the corresponding two- wave predictions (at £7/2, as shown by 
the dotted line), particularly at small k. 

Although it is apparent from Figure 0] that both \s c — s n \ and \s ac — s n \ vanish as U — > 0, we 
would really like to show that the relative error of the predictions vanishes as U — > 0, and hence that 
s c and s ac have the correct limiting behaviour at small amplitude. Therefore, given a numerically 
determined growth rate s n , for a theoretical prediction s p we introduce a relative error 

e r (s P ) = (61) 

The quantity |e r (s c )| is shown in Figures 2fc,d for the runs of Figures 2t,b. In both the two-wave 
and four- wave cases one can see that the relative error does vanish as U — > 0. Indeed, one can show 
in each case that e r w U 2 x some function of (k/U). Note that the relative errors for the two- wave 
and four-wave cases are of the same order of magnitude when comparing by U (rather than by the 
individual wave amplitude sCq). 
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Figure 4: Numerically determined growth rates and theoretical predictions at c = 5 = 1.5, k = 1 as a function 
of disturbance wavenumber kd = k . (a) Two- wave results, (b) Four- wave results, (c) The quantity |e r (s c )| 
for the two-wave prediction s c , as defined in JBIJ. (d) The quantity |e r (s c )| for the four-wave prediction s c . 
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For the two- wave cases we have already noted that s c and s ac are almost indistinguishable. How- 
ever, is s ac more accurate than s c ? Perhaps, for the larger amplitude cases in Figure EK, 63 is a more 
suitable model than (|33j) because the wave amplitudes are now of the same order as (v 91 —v 92 ) = 0.48? 
Or perhaps the coupled NLS model (|36|) is necessarily more accurate than the non-locally coupled 
NLS model (ffify? At least for these parameters, the answer to all three of these questions is no, since 
\s n s c \ is marginally smaller than \s n s ac |. Further, from (ffify we calculate that \v 92 — v 91 \ > (3 for 
all the two- wave cases, so that at the wavenumber k C2 where the uncoupled C-mode becomes unstable 
the coupled NLS equations predict the destabilizing scenario 1(a) of Section POl However, one can 
clearly see in the large amplitude cases of Figure EJi that there is a stabilization at k = k C2 , rather 
than a destabilization. 



The second case we discuss, also a full exchange solution, has c = VL5, 8 = 3, k = 4, for which 
v 91 = 1.07, v 92 = 0.95, so that v 91 — v 92 = 0.12. The nonlinear interaction coefficients are \i\ — 0.05, 
ji2 = 0.38, /i3 = 0.67, /i4 = 0.34, and uj'{ = —0.12, uj 2 = 0.19. The results are shown in Figure [5] For 
the two- wave case, both s c and s ac offer a good approximation to s n . At small amplitudes, s c and 
s ac are almost identical, as expected, and as shown in Figure Efc, s c shows the desired convergence 
properties as U — > 0. At larger amplitudes, where perhaps we enter the regime with U ~ (v 91 —v 92 \ 
there are differences between s c and s ac , and here s ac is a better approximant. Indeed, \v 92 — v 91 \ > (3 
as defined in so that the coupled NLS models predict the destabilizing scenario 1(a) of Section IOI 
with s n > at k = fc C2 , which is verified at ?7 = 0.6, for instance. For the four- wave case, there are 
clear differences between s c and s ac even at small amplitudes when k is also small. To obtain relative 
errors comparable to the two- wave case, it is now necessary to use the prediction s ac , as shown in 
Figure Eli. In both two- wave and four- wave cases, note that \e r \ increases with k at fixed U, which is 
perhaps associated with a breakdown in the scale assumption implicit in the weakly nonlinear theory 
leading to s c and s ac . 

Differences between the theoretical predictions s c and s ac can be highlighted by considering pa- 
rameters such that v 91 — v 92 . Therefore we consider the case c = 2/3, S = 3/4, k = 3/2, for 
which v 91 = v 92 = 0.71. The nonlinear interaction coefficients are fi\ = 0.07, ji2 = 0.29, n% = 0.60, 
/i4 = 0.57, and uj" = —0.18, (jJ 2 = 0.37. These do not satisfy the full exchange condition (fTTj) . For 
two-wave solutions, (|39|) gives sAo = 0.13/7, eCo = 0.37/7, and for four-wave solutions gives 
eAq = sBq = 0.065 Z7, eCo = sDq = 0.18 U. According to Q and (|T5j). the maximum amplitude of w 
in the linear solution is 1.17/7. 

Results for these parameters are shown in Figure El For both two- wave and four- wave cases, there 
are now significant differences between s c and s ac . Although s c gives moderate agreement with s n , 
it completely misses the second tongue of instability at large k, and |e r (s c )| does not decrease with 
U. However, s ac does have the correct qualitative behaviour. To assess its quantitative performance, 
we consider the relative error e r (s ac ), as shown in Figures EJ,d. Because s ac can be greater or less 
than s n , we plot e r (s ac ) (which may be positive or negative) rather than |e r (s ac )| as before, allowing 
the form of the curves to be fully appreciated. For the two- wave case, the points with k/U < 5/6, 
which includes the region of maximum growth, show good convergence as U — > 0. For k/U > 5/6, 
convergence is no longer clear, although the actual behaviour could be masked by the limited discrete 
resolution of e r . Nevertheless, errors here remain at about 10%. For the four- wave case, the points 
with k/U < 0.32, which includes the region of maximum growth, also show good convergence as 
U — > 0. For k/U > 0.32, it appears that e r (s ac ) assumes a form independent of Z7, with no decrease 
as U 0, although errors here still remain at about 10%. 

For the two- wave cases we have (3 > \v 92 — v 91 1 = as defined in (|l8j) . so that we lie in scenario 
1(b) of Section [Ql One can then verify that R < R c for all the cases in Figure Et? so that the 
coupled NLS models predict a destabilization at k — fc C2 , as is clearly the case at each U. Note this 
destabilization at k = k C2 is not inconsistent with the vanishing of s ac (and s n ) for some k < fc C2 , as 
noted at the end of scenario 1(b). This twin-tongued growth rate structure is typical for two- wave 
cases with \v 91 — v 92 \ <C 1. 

For the four- wave cases, the numerical simulations also show the existence of a second tongue of 
instability with a higher disturbance wavenumber than those shown in Figure E3 This appears at 
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Figure 5: Numerically determined growth rates and theoretical predictions at c = a/L5, 5 = 3, /c = 4 as a 
function of disturbance wavenumber kd = k. (a) Two-wave results, (b) Four-wave results, (c) The quantity 
\e r (s c ) | for the two- wave prediction s c , as defined in (|61|> . (d) The quantity \e r (s ac ) | for the four- wave prediction 
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U = 0.3 as an instability with s n = 0.00256 and k = 15/34, and at U = 0.4 as an instability with 
s n = 0.00483 and k = 3/7, and in both cases is the strongest instability observed at that amplitude U. 
Such instabilities are not predicted by any of our theories. Indeed, even though this is an instability of 
a weakly nonlinear plane wave, the modulating wavenumber k ~ fc/3, so that it does not correspond 
to a long spatial modulation of the type considered here. 

6.4 Character of the instability at finite amplitude 

In our previous study of the instability of two- wave solutions ^Hj ? we showed how the energy exchange 
between u and w in the plane-wave solution can have a counterpart during the finite amplitude stage 
of a modulational instability. We now examine further the nature of this finite amplitude energy 
exchange, and extend the discussion to the four-wave case and to the long-time behaviour of the 
system. 

In order to make clear conclusions, we consider a restricted class of scenarios. First, we restrict 
the plane- wave parameters to c = 5 and k = 1, all of which satisfy the full exchange condition 
Then, for c and S of 0(1), v 91 — v 92 w 0.5, and we lie in a stability regime like that of Figure SI 
Second, we restrict the computational geometry so that just one of the numerically resolved sidebands 
is unstable. We achieve this by taking I « 27r//c c *, where k c * is the most unstable wavenumber 
associated with the prediction s c , as determined by ((44)) or by the C-mode version of We then 
obtain an evolution in which the initial plane-wave structure becomes modulated on the lengthscale 
of the computational domain, with no possibility of an interaction between modulational instabilities 
of two different wavelengths. 
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Figure 6: Numerically determined growth rates and theoretical predictions at c = 2/3, S = 3/4, = 3/2 
as a function of disturbance wavenumber kd = k . (a) Two-wave results, (b) Four-wave results, (c) The 
quantity e r (s ac ) for the two- wave prediction s ac , as defined in (|61f) . (d) The quantity e r (s ac ) for the four- wave 
prediction s ac - 



We start by discussing an important special case, that of c = S = 1. Here, the system © 
decouples into a Sine-Gordon equation and a linear wave equation, and, therefore, is integrable |25| . 
As a consequence, the modulational instability of plane- wave solutions shows a recurrence, similar 
to the famous Fermi-Pasta-Ulam recurrence 41 j. We view this behaviour by considering the time 
evolution of three quantities: E p i anei as defined in (|59j) . representing the strength of the plane- wave 
and its nonlinear harmonics; E side \i, as defined in ((5%)) . representing the strength of the modulational 
instability in the first sideband pair, with disturbance wavenumber k = 2k /I; and E — E p i ane — E side \i, 
representing the contribution from other wavenumbers. 

These quantities are shown in the top row of Figure [3 during the instability of a two- wave solution 
at c = 5 = 1, with U = 0.4 and I — Sir. One can clearly see an almost periodic recurrence, with most of 
the energy either in E p i ane , or E side \\. In the second row, one can see how E is partitioned between u, 
w and the coupling, as defined in whilst in the third and four rows one can see u(x, t) and w(x, t) 
respectively (the complete computational domain is shown). At first the plane- wave solution has a 
periodic energy exchange, which is almost complete between E u and E w , as approximately described 
by H3 and (0- At the first minimum in E p i anei near t — 1300, there is still a periodic energy 
exchange between u and w, although the exchange is no longer almost complete, i.e. E u and E w no 
longer almost vanish periodically. Further, at this time u{x,t) and w{x,t) are strongly modulated 
on the largest resolved scale, and the large amplitude behaviour pulses between u and w. Thus, the 
almost complete energy exchange with spatially uniform structures initially observed is modified in 
the nonlinear regime to an incomplete energy exchange with structures localized on the lengthscale of 
the computational domain. Note that these same two scenarios alternate, in line with the recurrence 
of E p i ane and E si d e \i- One can see how u(x,i) and w(x,t) return almost to their initial plane- wave 
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Figure 7: Time evolution of a two-wave solution at c = 8 = 1, k = 1, U = 0.4, with / = Sty. 

form when E side \i is minimised, near t = 13, 170 for instance, whilst the pulsing modulations return 
when E side \i is maximised, near t = 17, 150 for instance. 

The four- wave case at c = S = 1 with U = 0.4 is similar, as shown in Figure |H1 Once again there 
is an almost periodic recurrence in E p i ane and E side \i. However, the partitioning of E between the 
components E Ul E w and E c is more complicated than in the two- wave case, with rapid oscillations 
in each component superimposed on a slower exchange between E u and E w . Once again, when 
Epiane is maximised the slow exchange between E u and E w is almost complete, whilst when E side \i 
is maximised the slow exchange becomes incomplete. The general pattern of an almost periodic 
sequence of alternating exchange patterns remains, from approximately uniform wavetrains to localized 
structures and back, and so on. 

A natural question to ask is to what extent similar behaviour occurs when 1. Shown in Figure[2l 
are results of corresponding simulations at c = S = 1.01 and c = 5 = 1.1, with U = 0.4 in both cases. 
For the two- wave cases, shown in panels (a) and (c), the period of the recurrence now varies, although 
the character of each recurrence event appears to remain almost the same. For the four-wave cases, 
shown in panels (b) and (d), the character of the recurrence events now varies, but the evolution is 
still predominantly characterized by an exchange between E p i ane and E side \i. Thus both two- wave 
and four- wave cases exhibit an imperfect recurrence, which we might expect since, for c close to 1, 
the equations are close to being an integrable system. Interestingly, the system (0 with c close to 1 
has been used to model some dynamical processes in the DNA [23J [22j ? wnere such energy exchange 
processes could, potentially, play an important role. 
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Figure 8: Time evolution of a four-wave solution at c = S = 1, k = 1, U = 0.4, with / = 16 7r. 
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Figure 9: Time evolution of solutions with k — 1 and U — 0.4: (a) a two- wave solution at c — 5 = 1.01, with 
/ = 8tt; (b) a four-wave solution at c = S = 1.01, with / = 167r; (c) a two-wave solution at c = S = 1.1 with 
/ = 8tt; (d) a four-wave solution at c = S = 1.1 with / = 16 7r. 
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Figure 10: Time evolution of a two-wave solution at c = 5 = 1.2, k = 1 and U = 0.4, with / = 10 7r. 

Although the recurrence becomes imperfect when c ^ 1, the alternating pattern of approximately 
uniform wavetrains and localized structures remains. Shown in Figures 1X771 and ITT1 are the results of 
two-wave and four-wave simulations at c = S = 1.2, with U = 0.4. The two-wave case of Figure ITUl 
still looks remarkably similar to that of Figure [3 taken at c = 5 = 1. However, in the four- wave case 
of Figure fTTl the recurrence breaks down around t = 13,000, only to reappear later in the evolution. 
Further, during the second sequence of recurrence events, around t = 17, 000, the alternating pattern of 
approximately uniform wavetrains and localized structures reappears. For instance, around t = 17, 100 
one can see that the solution only has a small spatial modulation, whilst around t = 17, 730 the solution 
has a spatial character almost identical to that of the first recurrence event around t = 5680. More 
remarkably, note from the second row of Figure ITT1 how even the details of the exchange between E u , 
E w and E c are almost identical around t = 5680 and t = 17, 730. 

However, for larger values of c = 5, evolutions which might be described as recurrent become rare, 
and typically the solutions display a more irregular behaviour. To illustrate this, we show results for 
the case c = 5 = 1.5. Shown in Figure IT2l is the evolution of a two-wave solution, with U = 0.4 and 
I = 12 7r, whilst shown in Figure IT31 is the evolution a four- wave solution, with U = 0.4 and I = 24 7r. 
(One can pick out these parameters in Figure and verify that they correspond to relatively strong 
instabilities in each case). In the two- wave case, at large times the solution loses much of the character 
of the earlier recurrence. The four- wave case displays an irregular, imperfect recurrence, the character 
of which does seem to be similar to the earlier events. 
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Figure 11: Time evolution of a four-wave solution at c = S = 1.2, k = 1 and U = 0.4, with I — 20 7r. 
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7 Conclusion 

We have studied the nonlinear dynamics of two pairs of counter-propagating waves within the frame- 
work of a two-component system of coupled Sine-Gordon equations (|3J). Using both analytical and 
numerical methods, we considered the evolution of spatially periodic solutions characterized by a sin- 
gle dominant wavenumber /c, focussing on the possible modulational instability of spatially uniform 
plane- wave solutions. We concentrated on a particular class of plane- wave solutions corresponding to 
an energy exchange between the two physical components of the system, and for these solutions we 
noted how modulational instabilities can lead to a modification of the energy exchange, and to the 
formation of localized structures. 

The analytical results were based upon asymptotic multiple-scales expansions for the evolution 
of weakly nonlinear wave packets, which reduced © to coupled equations for the evolution of four 
spatially varying wave amplitudes. However, the form of the amplitude equations depends crucially 
upon the assumed length-scale of the wave packets vis-a-vis the wave amplitude parameter e « 1 
and the group velocities v gi and v 92 . In Section [3| we derived two classes of evolution equations 
for the regime relevant to this system with v gi and v 92 both of 0(1): a set PT|) of non-dispersive 
nonlinearly coupled envelope equations for modulations on a lengthscale ~ £ -2 , and two sets and 
of non- locally coupled NLS equations for modulations on a lengthscale ~ £ _1 , according as to 
whether (v 91 — v 92 ) = 0(1) or 0(e). We have also considered a set (|55j> of coupled NLS equations 
with asymptotically weak dispersion, which formally includes each of the other models as different 
asymptotic reductions. However, it is not clear which set of amplitude equations most efficiently 
captures the nonlinear dynamics of interest. 

The stability of plane-wave solutions involving only a single pair of waves was considered in Sec- 
tion 21 whilst the stability of plane- wave solutions involving two pairs of counter-propagating waves 
was considered in Section In both cases the non-locally coupled NLS equations for the regime 
with (v 91 — v 92 ) = O(l) lead to quadratic dispersion relations for the disturbance frequency, enabling 
simple stability conclusions to be reached. One interesting feature of these results is that a coupled 
four- wave solution can be stable even whilst one or more of the four component waves can be unstable 
in isolation. This is in contrast to coupled two-wave solutions, for which stability is equivalent to that 
of the two component waves in isolation. For the regime with (v 91 — v 92 ) = 0(e), the non-locally cou- 
pled NLS equations lead to quartic dispersion relations for the disturbance frequency. The coupled 
NLS equations lead to a quartic dispersion relation for two-wave solutions, and to an eighth-order 
dispersion relation for four-wave solutions, although the latter was not evaluated here. 

Numerical methods were used to solve the fully nonlinear Sine-Gordon equations © in a periodic 
domain, using plane- wave initial conditions, as described in Section A major issue was to assess 
the performance of the stability predictions obtained from the different weakly nonlinear models. For 
cases with (v 91 — v 92 ) — 0(1), the quadratic dispersion relations of the non-locally coupled NLS 
equations give good predictions: relative errors are about 2 % when the amplitude parameter U = 0.5, 
and these decrease with U. For two-wave cases, the quartic dispersion relation of the coupled NLS 
equations gives no more accurate predictions, even at large amplitudes. As \v 91 — v 92 \ decreases, the 
quadratic dispersion relations become less accurate, and the quartic dispersion relations of the non- 
locally coupled NLS regime with (v 91 — v 92 ) = 0(e), or equivalently of the coupled NLS equations in 
the two-wave case, become more useful. For cases where the group velocities are very close or equal, 
the quadratic dispersion relations of the non-locally coupled NLS equations no longer even predict 
qualitatively the correct form of the dispersion curve, and a quartic dispersion relation must be used. 
The region of maximum growth is then well described, with relative errors decreasing with U, although 
at larger disturbance wavenumbers convergence appears to be poor and relative errors are about 10 %. 
In some four-wave cases there is an additional dominant short-wavelength instability, which does not 
correspond to a long spatial modulation of the type considered theoretically, and is not predicted by 
any of our theories. 

The numerical simulations also enabled us to examine the long-time evolution of modulational 
instabilities. For the energy exchange parameters of interest, the initial solutions are spatially uniform 
plane- waves with an almost complete energy exchange between the components, whilst later in the 
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evolution the solutions can become spatially localized with an incomplete energy exchange. Further, 
when the system is close to the integrable case at c — 1, the time evolution is distinguished by a 
remarkable almost periodic sequence of energy exchange scenarios, with spatial patterns alternating 
between approximately uniform wavetrains and localized structures. As c moves away from 1, the 
time evolution becomes more disordered. 

Although the whole of this study was performed for a two-component set of coupled Sine-Gordon 
equations, we would expect similar phenomena to arise in a wide class of nonlinear wave systems. Thus, 
in addition to applications to multi-component systems such as those mentioned in the Introduction, 
processes of the type considered here should arise in other physical systems which can support several 
wave modes. For instance, an interesting example might be energy exchange between nonlinear internal 
waves in a layered, or continuously stratified, fluid. 
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